A fast algorithm for the hnear canonical transform 



Rafael G. Campos and Jared Figueroa 
Facultad de Ciencias Fisico-Matematicas, 
Universidad Michoacana, 
58060, Morelia, Mich., Mexico. 

rcampos@umich.mx, jared@fismat.umich.mx 



MSG: 65T50; 44A15; 65D32 

Keywords: Linear Canonical Transform, Fractional Fourier Transform, Quadrature, Her- 
mite polynomials. Fractional Discrete Fourier Transform, f ft 



Abstract 

In recent years there has been a renewed interest in finding fast algorithms to compute 
accurately the linear canonical transform (LCT) of a given function. This is driven by the 
large number of applications of the LCT in optics and signal processing. The well-known 
integral transforms: Fourier, fractional Fourier, bilateral Laplace and Fresnel transforms 
are special cases of the LCT. In this paper we obtain an 0{N log N) algorithm to compute 
the LCT by using a chirp-FFT-chirp transformation yielded by a convergent quadrature 
formula for the fractional Fourier transform. This formula gives a unitary discrete LCT 
in closed form. In the case of the fractional Fourier transform the algorithm computes 
this transform for arbitrary complex values inside the unitary circle and not only at the 
boundary. In the case of the ordinary Fourier transform the algorithm improves the 
output of the FFT. 
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1 Introduction 



The Linear Canonical Transform (LCT) of a given function f{x) is a three-parameter 
integral transform that was obtained in connection with canonical transformations in 
Quantum Mechanics [H [2] • It is defined by 

\/2nib J-oo 

for b ^ 0, and by ^fde^'^y^ f{dy), if 6 = 0. The four parameters a, 6, c and d appearing in 
([1]), are the elements of a 2 x 2 matrix with unit determinant, i.e., ad — bc= 1. Therefore, 
only three parameters are free. Since this transform is a useful tool for signal processing 
and optical analysis, its study and direct computation in digital computers has become 
an important issue [3]- [10], particularly, fast algorithms to compute the linear canonical 
transform have been devised [H [7] . These algorithms use the following related ideas: (a) 
use of the periodicity and shifting properties of the discrete LCT to break down the orig- 
inal matrix into smaller matrices as the FFT does with the DFT, (b) decomposition of 
the LCT into a chirp-FFT-scaling transformation and (c) decomposition of the LCT into 
a fractional Fourier transform followed by a scaling-chirp multiplication. All of these are 
algorithms of 0(iVlogA^) complexity. 

In this paper we present an algorithm that takes 0{N log N) time based in the decomposi- 
tion of the LCT into a scaling-chirp-DFT-chirp-scaling transformation, obtained by using 
a quadrature formula of the continuous Fourier transform pTl IT2] . Here, DFT stands 
for the standard discrete Fourier transform. To distinguish this discretization from other 
implementations, we call it the extended Fourier Transform (XFT). Thus, the quadrature 
from which the XFT is obtained, uses some asymptotic properties of the Hermite polyno- 
mials and yields a fast algorithm to compute the Fourier transform, the fractional Fourier 
transform and therefore, the LCT. The quadrature formula is (9 (1/A^) -convergent to the 
continuous Fourier transform for certain class of functions [I3j . 



2 A discrete fractional Fourier transform 

In previous work [12], [13], [H], we derived a quadrature formula for the continuous 
Fourier transform which yields an accurate discrete Fourier transform. For the sake of 
completeness we give in this section a brief review of the main steps to obtain this formula. 
Let us consider the family of Hermite polynomials Hn{x), n = 0,1, . . ., which satisfies the 
recurrence equation 

H^+iix) + 2nH^_i{x) = 2xHn{x), (1) 
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with = 0. Note that the recurrence equation ([T]) can be written as the eigenvalue 

problem 

/O 1/2 •■■\ /Ho{x)\ /Hoix)\ 



1 1/2 
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Let us now consider the eigenproblem associated to the principal submatrix of dimension 
ATof (H 
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It is convenient to symmetrize Ti by using the similarity transformation ST-CS~^ where S 
is the diagonal matrix 

S = diag < 1, —;=, . . . , — , 

\ V2 ^{N-l)\2^-^ 

Thus, the symmetric matrix H = STiS^^ takes the form 

/o 
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The recurrence equation ([T]) and the Christoffel-Darboux formula [15] can be used to solve 
the eigenproblem 

Huk = XkUk, k = 1,2, . . . ,N, 

which is a finite-dimensional version of ([2]). The eigenvalues Xk are the zeros of H]\f{x) 
and the fcth eigenvector Uk is given by 

where Si, . . . , s^r are the diagonal elements of S and is a normalization constant that 
can be determined from the condition u]^ = 1, i.e., from 

2 Hn{Xk)Hn{Xk) _ ^ 

2^ - ■■ - J-- 



n=0 



2"n! 
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Therefore, 

'2^-1 (A^-l)! (-1)^^+^ 



Ck 



Thus, the components of the orthonormal vectors Uk, k = 1,2, N, are 



y N{n- 1)1 HN-i{xk) 

n = 1, . . . ,N. Let U be the orthogonal matrix whose fcth column is Uk and let us define 
the matrix 

J^, = V2^U-^D{z)U, 

where D(z) is the diagonal matrix D{z) = diag{l, 2;, 2;^, . . . ,2:^"^} and z is an complex 
number. Therefore, the components of J-'z are given by 

( —Di+k oN-i ( ivf _ T\] 7" 
N HN-.i{xj)HN-i{xk) ^ 2"ra! 

Next, we want to prove that if is large enough, (jl]) approaches the kernel of the fractional 
Fourier transform evaluated ai x = xj, y = Xk- To this, we use the asymptotic expression 
for HNix) L15J) 

^^^^^ " r(AV2 + l) ^°<v'2iVTT X - ^). (5) 
Thus, the asymptotic form of the zeros of Hi\f{x) are 

/2A;- A^- 1\ vr 

= [ V2N ) r 

A; = 1, 2, . . . , iV. The use of (jS]) and §^ yields 



00, 



and the substitution of this asymptotic expression in (jl]) yields 

^ n=0 

Finally, Stirling's formula and Mehler's formula pjj produce 



1 2 (l + z^)(^^+^|)-4a:j;>;feZ 

^ y e ^^^^^^ Axfe, (7) 

where Axk is the difference between two consecutive asymptotic Hermite zeros, i.e., 

A.» = .,„-x. = -|=. (8) 

4 



Let us consider now the vector of samples of a given function f{x) 

/=(/(Xi),/(x2),...,/(x^)r. 

The multiphcation of the matrix JF^ by the vector / gives the vector g with entries 

k=l * k=l 

where j = 1,2 ... ,N. This equation is a Riemann sum for the integral 



°° (l + z^)(y^+x^)-4xyz 



^z[f{x),y] = \l- I e f{x)dx, 



where \z\ < 1. Therefore, if we make yj = Xj, 



oo 



N 



k=l 



Note that ^z[g{x),y] is the continuous fractional Fourier transform [17] of g{x) except 
for a constant and therefore, J-'z is a discrete fractional Fourier transform. 



3 A fast linear canonical transform 

Firstly, note that if 6 7^ 0, the LCT can be written as a chirp-FT-chirp transform 



fOO 



-00 



Thus, for b 0, the LCT of the function f{x) can be represented by the 1/6-scaled Fourier 

2 idiP' 

transform of the function gix) = e^fc~/(x), multiplied by 



/2nib ' 



On the other hand, note that for the case z = ±i, ([7]) yields a discrete Fourier transform 
{T±i)jk — e^**-'*'' Atfc, that can be related to the standard DFT as follows. The use of ([6]) 
yields 

(^^)^., = e±^*^*^A4 = -^e'^^^-"^)^^-"^) (10) 

where we have used (j6]) and ([8]). Since Yl!k=ii.-^i)ikf{xk) is a quadrature and therefore, 
an approximation of 

/oo 
e'y^''f{x)dx, 
-00 

a scaled Fourier transform 

/•oo 

e'^y^^f{x)dt = giKy,) (11) 



has the quadrature '^k=i^jkfi^k), where 



Fjk 



(12) 



If we choose k, = A/tt, (|T2l) takes the form 



Fjk 



(JV-I)^ 

vre 2 AT 



e N " 



for j,k = 0,1,2, . . . , N — 1, and J2k=i^jkfi^k) is an approximation of g^Ayj/ir). If 
now we choose k, = Ab/n, but we keep the same matrix (fT3|l . then Yl!k=i^ikf{xk) is an 
approximation of 



If now we replace /(x) by e'aii /(x) and take into account ffTOl) . we have that 



AT 



A;=l 



is an approximation of the product of functions 



/2iTib 



-1 



£{a,b,c,4[j(^^)^^] evaluated at 



i/j = Abxj/iT. Therefore, a discrete (scaled) linear canonical transform L can be given in 
closed form. If we denote by G{y) the LCT of f{x), then 



N 



G{y,) = G{4bx,/7r) = J2{SiFS2)jkf{xk), 



k=l 



where 5*1 and 5*2 are diagonal matrices whose diagonal elements are e j ^^'Kib, and 
respectively. As it can be seen, the matrix L = S1FS2, which gives the discrete 



iaxj/2b 



e 1 



LTC, i.e., the XFT, consists in a chirp-DFT-chirp transformation, where DFT stands for 
the standard discrete Fourier transform. Therefore, we can use any FFT to give a fast 
computation of the linear canonical transform G = Lf. 

Now, the fast algorithm for the linear canonical transform can be given straightforwardly. 
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Algorithm 



To compute an approximation G = {Gi, G2, ■ ■ ■ , Gn)'^ of the hnear canonical trans- 
form G{y) = C^"-'^'^''^^[f{x),y] evaluated at yj = Abxj/n, where xj = (^^^^^=^^ |. 

1. For given set up the vector u of components 



-i^LtzlMjzll iax^ /2b r f 2k - N - 1 

Uk = e N e f vr ^= — 

V 2V2N 



k=l,2,...,N. 



2. Set yj = Abxj/iT and compute the diagonal matrix S according to 



Tie 2 N C 2b C N " ' 



j,k = 1,2,..., N. 



3. Let Dp be the discrete Fourier transform, i.e., {Dp)jk = e^^^^, j,k = 
0,1,2, N ~ 1. Obtain the approximation Gj to G{^Xj) by computing 
the matrix-vector product 

G = SDpu, (13) 

with a standard FFT algorithm. 



4 Example 

For this example we take an integral formula found in [18] that gives 



Giy) = , , e^^et-^-(^)e- ""4^::^ ' (14) 



j aey^ 2b(a^dy^+2f3ay + l3'^a) 

X e 2(462^2+^2) g 4(,2a2+a2 ^ 

if f[x) = e-("f'+2/3s/+7)^ a > 0. Figure 1 shows the exact LTC with a = 1, /? = 2, 7 = 3, 
a = l,6 = 2,c = l/2, and d = 2, compared with the approximation given by the XFT. 
Figure 2 shows the exact Fresnel transform of f{x) = e"*-"^ +2/3j/+7) -v^^ith a = 2, P = 1, 
7 = 3, a = 1, 6 = 100, c = 0, and d = 1, compared with the approximation given by the 
XFT. 



7 




-6 -4 -2 2 4 6 -6 -4 -2 2 4 6 

y y 



Figure 1: Real part (A) and imaginary part (B) of the exact linear canonical transform 
(solid line) compared with the output of the XFT (dashed line) computed with = 512. 
The function Giy) is that given in f|T4l) for a = 1, /5 = 2, 7 = 3, a = 1, 5 = 2, c = 1/2, 
and d = 2. 




Figure 2: Real part (A) and imaginary part (B) of the exact Fresnel transform (solid 
line) compared with the output of the XFT (dashed line) computed with = 1024. The 
function Giy) is that given in for a = 2, /3 = 1, 7 = 3, a = 1, 6 = 100, c = 0, and 
d=l. 
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5 Conclusion 



We have obtained a discrete linear canonical transform and a fast algorithm to compute 
this transform by projecting the space of functions onto a vector space spanned by a 
finite number of Hermite functions. The XFT is a discrete LCT given by a unitary 
matrix in a closed form in which the DFT can be found at the core, surrounded by 
diagonal transformations, which makes easy to implement it in a fast algorithm. Since 
this discrete LCT is related to a quadrature formula of the fractional Fourier transform, 
it yields accurate results. 
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